NASA Contractor Report 3526 


A Variational Principle for 
Compressible Fluid Mechanics 

Discussion of the One-Dimensional Theory 



Robert Joel Prozan 


i o 


CC * > v'; fyS* TO 

; : L "’ 1 :■ 

S. * 4. 


CONTRACT NAS 1-16646 
APRIL 1982 





TECH LIBRARY KAFB, NM 



TECH LIBRARY KAFB, NM 



0Db2214 


NASA Contractor Report 3526 


A Variational Principle for 
Compressible Fluid Mechanics 

Discussion of the One-Dimensional Theory 


Robert Joel Prozan 

Continuum, Inc. 
Huntsville, Alabama 


Prepared for 

Langley Research Center 

under Contract NAS 1-16646 


f\JASA 

National Aeronautics 
and Space Administration 

Scientific and Technical 
Information Branch 


1982 



1.0 SUMMARY AND INTRODUCTION 


Conventional numerical analyses of fluid dynamic systems 
variously use finite element or finite difference analogs 
to satisfy the governing equations of motion. To date, 
numerical solutions are more of an art than a science. 

In this report it is proposed that the second law of 
thermodynamics, a neglected principle of physics, is of 
paramount importance to a successful solution to the 
governing equations of motion. The second law is not new, 
nor is the post priori use of the second law in ascertain- 
ing the physical correctness of the solution. What is 
new is the incorporation of the second law as a fundamental 
statement of motion which is of equal importance to the 
solution as are the familiar conservation equations. It 
is hypothesized that the second law of thermodynamics is 
actually a variational principle and, in point of fact, is 
used as such in the subsequent discussion. 

Many years of experience with the frustrations of numerical 
modeling have convinced the author that something is miss- 
ing in the conventional numerical approaches. Even if the 
conservation equations are satisfied the solution may fail. 
There must be other governing criterion that have hereto- 
fore been lacking. 

The ensuing discussion will use the second law of thermo- 
dynamics as a variational statement to derive a numerical 
procedure which provides insight to some fundamental 
questions not previously resolved. The procedure, based 
on numerical experimentation, appears to be stable pro- 
vided the CFL condition is satisfied. This stability is 
manifested no matter how severe the gradients (compression 
or expansion) are in the flow field. 

For reasons of simplicity only one-dimensional inviscid 
compressible unsteady flow will be discussed here; however, 
the concepts and techniques are not restricted to one 
dimension nor are they restricted to inviscid non-reacting 
flow. The solution here is explicit in time. Further 
study is required to determine the impact of the varia- 
tional principle on implicit algorithms. 
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2.0 DISCUSSION 


2.1 The Variational Principle 

In the subsequent discussion vector notation will be Used 
rather than the more frequently used Einstein notation. 
The equations of motion for compressible, inviscid, 
unsteady flow of a perfect gas are; 
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In text definitions will be used throughout. The symbols 
introduced above are: p (density) ; u (axial velocity); 

T (temperature) ; p(pressure); t(time); x(space coordinate); 
Z (gas constant); y (ratio of specific heats); E (specific 
internal + kinetic energy); H (total enthalpy). 

The Second Law of Thermodynamics is written; 
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The integral of the above equation over the region of 
interest and over an interval of time is; 

il[ / v ( IF + OF' dvdt - 0 

Equation (3) may be further integrated to yield; 
t 
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Thus the Second Law states that the entropy generated 
internal to the region under investigation must be station- 
ary or increase. Another view of the Second Law results 
when equation (3) is expanded; 
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Equation (5) states that if the conservation equations are 
satisfied exactly everywhere then the system entropy 
generation is zero. In a numerical analysis, however, the 
equations of motion are only satisfied approximately so 
that in general the equality will not be satisfied. The 
Second Law is therefore interpreted as stating that the 
"numerical" entropy created by the inexact scheme must 
be positive. There are many choices made during the con- 
struction of a numerical analog. Whether by accident or 
design, these choices must not (in the long run at least) 
violate the Second Law. 

The equations of motion which must be satisfied can be 
viewed as constraints to the entropy formation. Therefore, 
the constrained entropy function (S) can be written as; 

s = Jt^< / V ( ftT + X ‘ ( !f + !i> )dv + / cs PUs-dA)dt (6) 
where the A=(A , X , A „) are the Lagrange multipliers. 

P P U P £j 

Note that the constraint term looks very much like the 
familiar Method of Weighted Residuals (MWR) . Now let the 
time interval be small and let At=t 2 -t 1 . The one-dimen- 
sional domain will consist of N nodes which, for simplicity, 
will be subdivided into (N-l) elements of length Ax. 

(See figure 1 below. ) 
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Figure 1 Schematic of region 
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The numbering system is such that the element (e) is 
bounded by nodes (n) and (n+1) or that element (n) is 
bounded by nodes (n) and (n+1) . 

We will first develop the constraint term ($) . Let t=AtT 

and x-x =Ax£. The constraint term is then; 
n 


T=1 E 
$ = At Ax/ y 

t= 0 e=l 


A* (tt + 7— |f)dCdx 
J ^ = Q At 3 t Ax 3 E, 


(7) 


For the element shown in figure 2 below 
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Figure 2 Schematic of element 


the functions f and g are assumed to vary linearly. That is 

f = (1-OV^n+i' 9 = d-Ogn+^n+l* Then 
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We will also assume a linear functional for X viz; 
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(9) 


In the above expression the superscript ( ) is used to 
denote a dyad. Only diagonal dyads will occur during 
this development so that a vector pre or post multiplying 
the dyad will give the same result. There are also some 
occasions when it is necessary to redefine a vector as a 
dyad. It is felt that no undue interpretation problems 
to the reader should result. 

Now A ]_ n ~Ai n (£) and ^2n = ^2n^^* These functions are assoc- 
iated with the element (e) (left node (n) and right node (n+1) ) 
The constraint term becomes; 
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where (the time integration has been performed) , 


t n = fS +1 -fi ?; fn+1 = fgti-f^ +1 
At At 

and where the superscripts denote the time level 


Since the stated objective was to create an explicit 
scheme (this has nothing to do with the variational state- 
ment but is of great practical importance) , the weight 
functions (to use MWR terminology) are chosen such that; 

5=1 = S=l = 

/ A MS = 0; / A (l-5)d5 = 0. 

5=0 n 5=0 * 

If A, , A n are linear we find that; 
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where b-^ n , b 2 n are arbitrary constants. Now; 
^ =Q (A ln (1- ^ ; A 2n^ ; A ln' A 2n* = (-b ln ' b 2n ' " b ln 

— ■ I = — I = 
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Since the b were arbitrary constants so are the b. 
constraint term is now written as; 


The 
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n ± 2 n-1 n n n Ax n-1 Ax n 


where, in order £o achieve, a more recognizable form a 
substitution of b 2 n = (l-$ n ) ; b^ n = -$ n was made. 
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The variational statement of equation (6) becomes 


T = 1 

s = / / 

T = 0 V 
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(12) 


Differentiation of equation (12) with respect to the nodal 
Lagrange multipliers yields the equations of constraint; 
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for 2<n<N-l 


The above equations are an explicit numerical analog to the 
governing^ (or constraint) equations of motion. ^To deter- 
mine the IT , however, the arbitrary parameters S must be 
found. LeS, for example, 3 n =0 for all n. Then the interior 
nodal analog becomes; 


^n+ 2n S n=1 = 0 (14a) 


while if $ n =l we get; 

= 0 (14b) 


which are backward and forward differencing schemes 
respectively. The @ n therefore control the type = of 
analog generated. It is also apparent that the 3^ 
control the portion of each element that is associated 
with the node. Figure 3 illustrates the point. 


(13a) 

(13b) 

(13c) 
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Figure 3 Effective control volume 


The control volume is seen to be Av =A x (1-g ^ +3 n ) . 

Returning to the constrained equation (12) ana performing 
the integration on the same control volume basis yields; 
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where the flux term was expanded first order in time and 
integrated with respect to time. Differentiation of the 
above functional with respect to the T at each node yields 
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Equations (16) therefore define the Lagrange multipliers 
so that, for an interior node; 
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(16c) 


( 17 ) 
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Finally, using the above results, the constrained functional 
is differentiated with respect to the 3 n - The result is; 


3S 

33 
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= ixAt ( (|£Sfl±I_|£S a) . (g n+1 -g„) ) 


3f 


(18) 


Equation (18) indicates that for this low order analysis 
it is not possible to find a stationary value of the 
functional. This is of course due to the fact that the 
functional is linear in 3 n » What can be done is to determine 
the value of 3 n which is in the direction of increasing 
system entropy. If 3n is viewed as an interpolant then it 
must have a value between zero and unity. The approach 
taken for the purposes of the following computations is to 
choose 3 n equal to unity if the term in equations (18) is 
positive and equal to zero if the term is negative. In this 
fashion the differences become one-sided in the direction 
of maximum increase (or minimum decrease) of system entropy. 
It is also interesting to note that the boundary conditions 
naturally evolve from the analysis. Suppose that the left- 
most node remains at a fixed condition for all variables. 

In this case 6^ is not a variable and may not be determined 
by equation (18) , rather it must be zero so that (13a) is 
removed from the equation list. In the diaphragm burst 
problem to be discussed later, the axial velocity must be 
fixed (zero) while density and energy are variable. 

Therefore the value of 3 n for the momentum equation must be 
zero while the other two are determined by equation (18) . 

It can be seen, therefore, that the analysis determines a 
differencing scheme for each equation at each node for each 
time step which simultaneously satisfies boundary conditions 
and the condition of maximum entropy increase. As will be 
seen, depending on the values of 3n an Y particular equation 
at any particular node may be differenced forward or back- 
ward or centered or unchanged. 

Equations (13) may be solved for the f at each node. This 
results in a nonconservative analog. If, however, 
equations (13) are solved in the following fashion; 


Ax3]/fi = -8]_* (g2“9i) 

Ax (1-3 N _]_) • % = - (l~8 N _i) • (g N -g N _j_) 

Ax(l-3 n -i+l n ) ‘ f n = -I n ’ <5n+l~in> " (1 - = ?n-l ) * ^n“5n-l> 


(19a) 

(19b) 

(19c) 
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then the result is absolutely conservative. To illustrate 
this, note that; 


/ f dV = 
V 


N . 


I f n Av n = -<gN“9l> 
n-1 


( 20 ) 


so that all the interior flux terms cancel and just the 
entering and exiting flux terms remain. 

To summarize the computational procedure, having known 
values of f n at a given time level, then 

(1) compute gn = 9n^n) f° r nodes. 

(2) compute 9ps/9f for all equations at all nodes. 

(3) use equations X 18 ) (and boundary conditions) to 
determine the 3n f° r each element. 

(4) solve equations (19) for (f n AV n ) at all nodes. 

(5) integrate for next time level; 

~k+l _k , 

f n = f n + (fn Av n> At / Av 


(6) repeat (1) - (5) . 


2.2 Results 

The one-dimensional equations restrict the problems which 
can be solved to demonstrate the method. Two relatively 
difficult problems may be investigated, however. These 
are the shock tube and diaphragm burst. In the former 
case the boundary conditions dictate that Fi=( 0,0,0) while 
Fm_i are free. In the latter case Fi^F^-i 22 ( * r 0 , *) where 
the (*) indicates a free parameter. In all cases y=1.4 
and Z-l. 

Figure 4 illustrates the behavior of the method for the 
shock tube conditions. The initial disturbance of the 
square wave has not yet washed out. This impulsive start 
typically gives rise to wavelets which migrate back and 
forth between the inlet and the traveling primary wave. 

If the solution domain were long enough a better descrip- 
tion of the wave characteristics would result. In a higher 
spatial dimension analysis, however, economics prohibit 
the use of considerably more than the 40 points in a single 
direction. The traveling wave description shown in figure 4 
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is therefore typical of what would occur in a two-dimen- 
sional or three-dimensional analysis. The final velocity 
and pressure distributions are, of course, correct. 

Figure 5 illustrates typical results on a diaphragm burst 
problem. The results shown occur just before the shock 
wave impinges on the downstream wall and just before the 
expansion wave reaches the upstream wall. Calculations 
have been run for many thousands of time steps in which 
the waves reflect many times from both ends. These results 
are not particularly interesting and in the interest of 
brevity are not shown here. 

A typical time step is shown in figure 6. The rightmost 
three columns contain the values of beta for each equation 
at each node for this time step. These values of beta 
correspond to differencing schemes which are identified in 
the column headed by DIFF. The three letters at each node 
identify the scheme constructed for the three equations of 
motion. U means unchanged, B means backward, F means for- 
ward and C means centered. Time steps exceeding the CFL 
condition, as well as negative pressures and/or densities, 
are physically unrealistic. In the context of the varia- 
tional approach they must be viewed as inequality constraints. 
Due to the difficulty of formally entering an inequality 
constraint into the functional definition, these constraints 
were handled in the following fashion: The time step was 

never allowed to exceed some predetermined fraction of the 
least local CFL and should a violation of the pressure or 
density inequality constraint occur then the time step was 
halved. Subsequent time steps were increased by 10 percent 
each step until the maximum time step is reached. It seems 
reasonable to expect that, particularly during initial 
transients, time steps may be attempted which would lead 
to negative pressures or densities. This should only be a 
temporary condition. Instability would be manifested by 
repeated attempts to achieve unrealistic values; fortun- 
ately this never occurred. 
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SHOCK TUBE ANALYSIS 


Initial Conditions: p=l,p=l,u=0 (field); p-5 , p-2 . 5 ,u-2 (inlet) 
Boundary Conditions: fixed inlet (left end); free outlet 



Figure 4 Typical shock tube solution 


nressure 


DIAPHRAGM BURST ANALYSIS 



Initial Conditions: p=l , p=l ,u=0 ; x>20 p=5 , p=2 . 5 ,u=0 ; x<20 

Boundary Conditions: ends closed 



Figure 5 Typical diaphragm burst solution 
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Figure 6 Illustration of differencing 



3.0 CONCLUSIONS 


In the preceding discussion it was proposed that the 
Second Law of Thermodynamics is a variational principle 
applicable to compressible gas dynamics. It was shown 
that the principle states that the system entropy must 
not decrease and that the equations of motion act as 
constraints to the entropy formation. It was further 
shown that, at least in an explicit formulation, arbi- 
trary parameters exist which control or dictate the 
type of differencing scheme. The analysis then discusses 
how to select these arbitrary parameters such that the 
variational principle is satisfied while maintaining 
exact conservation. 

The concept is demonstrated numerically for the shock 
tube and diaphragm burst situations. The success of these 
calculations demonstrates the validity of the original 
proposition, i.e. , that the Second Law is indeed an 
applicable variational principle and furthermore demon- 
strates the validity of the calculational procedure 
derived from the stated principle. 

A fundamental advance in the state of the art of comp- 
utational fluid mechanics can be expected as a result of 
these findings. 
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